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The behavior of orbits in charged-particle beam transport systems, inchiding both hnear and 
circular accelerators as weh as final focus sections and spectrometers, can depend sensitively on 
nonlinear fringe-field and high-order-multipole effects in the various beam-hne elements. The inclu- 
sion of these effects requires a detailed and realistic model of the interior and fringe fields, including 
their high spatial derivatives. A collection of surface fitting methods has been developed for ex- 
tracting this information accurately from 3-dimensional field data on a grid, as provided by various 
3-dimensional finite-element field codes. Based on these realistic field models, Lie or other methods 
may be used to compute accurate design orbits and accurate transfer maps about these orbits. Part 
I of this work presents a treatment of straight-axis magnetic elements, while Part II will treat bend- 
ing dipoles with large sagitta. An exactly-soluble but numerically challenging model field is used to 
provide a rigorous collection of performance benchmarks. 



I. INTRODUCTION 

For the design of high-performance linear accelerators, 
synchrotrons, and storage and damping rings it is es- 
sential to have realistic electric and magnetic field infor- 
mation for the various beam-line elements in order to 
compute accurate design orbits and accurate high-order 
transfer maps about the design orbits. There are similar 
requirements for other charged-particle beam transport 
systems such as final focus sections, high-resolution spec- 
trometers, and high-resolution electron microscopes. 

Realistic field data can be provided on a grid with the 
aid of various 3-dimensional finite element codes, some- 
times spot checked against measured data. But the com- 
putation of high-order transfer maps based on this data 
appears to pose an insurmountable problem: the calcula- 
tion of high-order transfer maps requires a knowledge of 
high derivatives of the field data. The direct calculation 
of high derivatives based only on grid data is intolera- 
bly sensitive to noise (due to truncation or round-off) in 
the grid data H]. We will see that this problem can be 
solved by the use of surface methods. The effect of nu- 
merical noise can be overcome by fitting field data on a 
bounding surface far from the beam axis and continu- 
ing inward using the Maxwell equations. (We recall that 
interior fields are uniquely specified by their values on 
bounding surfaces.) While the process of differentiation 
serves to amplify the effect of numerical noise, the pro- 
cess of continuing inward using the Maxwell equations is 
smoothing. This smoothing is related to the fact that har- 
monic functions take their extrema on boundaries. When 
using surface methods, all fits are made on such bound- 
aries. Therefore, if these fits are accurate, interior data 
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based on these fits will be even more accurate. 

In this paper we will devote our attention to magnetic 
beam-line elements. Static electric beam-line elements, 
and static electric and magnetic beam-line elements such 
as velocity selectors (Wien filters), could be treated in a 
similar way. For a treatment of RF cavities, see 

There are two magnetic cases that it is convenient to 
handle separately: straight and curved. For straight 
magnetic elements such as solenoids, quadrupoles, sex- 
tupoles, octupoles, etc., and wigglers, it is convenient to 
employ cylindrical surfaces. For the case of curved mag- 
netic elements, such as dipoles with large design-orbit 
sagitta, it is convenient to employ the surface of a bent 
box with straight ends. In all cases the bounding sur- 
face will surround the design orbit within the beam-line 
element and will extend into the fringe-field regions out- 
side the beam-line element, thus taking into account all 
fringe-field effects as well as all effects within the body of 
the beam-line element. 

For the case of straight beam-line elements it is con- 
venient to describe the magnetic field in terms of a mag- 
netic scalar potential -0. Then, if one wishes to compute 
transfer maps in terms of canonical coordinates, one can 
proceed with the aid of an associated vector potential 
A computed from ■;/;. Alternatively, if one wishes to in- 
tegrate noncanonical equations employing the magnetic 
field B, it can be obtained from the relation B = Wip. 

For the case of curved beam-line elements it is conve- 
nient to work directly with the vector potential. Its use 
in the case of canonical coordinates is then immediate. 
If instead one wishes to integrate noncanonical equations 
employing the magnetic field B, it can be obtained from 
the relation B = V x A. 

In this paper we will treat the case of straight beam- 
line elements. For this purpose we will employ the sur- 
face of a (virtual) cylinder of uniform cross-section (cir- 
cular, elliptical, or rectangular) surrounding the beam, 
and which lies within all iron or other magnetic sources. 
In these Green function can be determined for 

the geometry of the fitted domain as a series composed 
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of known, orthogonal, and complete functions. In the 
case of bent-bore magnetic elements, such as dipoles with 
large sagitta, it is not possible to obtain a suitable Green 
function. New techniques have therefore been developed 
for treating cases of such general geometries, and these 
techniques will be presented in a subsequent paper. 

Our emphasis will be on the accurate representation 
of fields in terms of scalar and vector potentials. Once 
these representations are available, there are a variety of 
methods for computing the associated design orbits and 
transfer maps about these orbits. An Appendix summa- 
rizes how this can be done when canonical coordinates 
are employed and the associated Lie-algebraic structure 
is exploited. 

A brief outline of this paper is as follows: Section II 
reviews circular cylinder harmonic expansions and intro- 
duces a collection of functions, known as on- axis gra- 
dients, which uniquely characterize the magnetic scalar 
potential. It also describes how the vector potential can 
be obtained from the scalar potential. The on-axis gra- 
dients themselves are generally unspecified functions of 
z. In some simple cases they can be found analytically. 
However, in general they must be determined numeri- 
cally. Section III describes how this can be done us- 
ing known magnetic field values determined at points 
on some regular 3-dimensional grid. In Section IV, we 
treat an analytically soluble model problem, that of a 
magnetic monopole doublet, which is used to benchmark 
these methods. In Section V, we discuss the advantages 
of these surface-fitting methods that result from numer- 
ical smoothing. In Section VI, we use these methods to 
study a proposed ILC damping ring wiggler. The paper 
concludes with a summary and an Appendix. Further 
detail may be found in [31 S] . 



II. CYLINDRICAL HARMONIC EXPANSIONS 



p'" cosm<j> = ne[{x + iy)"'], 



p™ sinmcf) = Im[{x -i- iy)"^]. 



(4b) 



(4c) 



We see that even powers of p and the combinations 
p™ cos mcf) and p™ sin mi/) are analytic (in fact, polyno- 
mial) functions of x and y. 

A general solution -0 satisfying the Laplace equation 
and analytic near the axis p = takes the form 



poo 

i;{p, 4>,z)^Y. * Im{kp)e'''' [G^,s{k) 
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+Gm,c{k) COS m0] 
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Note that the term containing Go_s{k) does not con- 
tribute to the above sum, and we may set Go,s = with- 
out loss of generality. By utilizing the Taylor series of the 
modified Bessel function /„, we may write ip in the form 
of a circular cylinder harmonic (multipole) expansion: 



lp{p,(j),z)^ ^ ['0m,s(p, 2;)sinTO0-|-'(/'m,c(p, 2:)C0STO(/)] , 

(6) 



m=0 

where for a ~ s, c. 
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The functions Cm]a, known as on-axis gradients [31 3], 
are defined by 



2'"m! 



& 6**^^/0™+" G™,„(fc). (8) 



A. Scalar Potential 

In a current-free region the magnetic field B is curl 
free, and can therefore be described most simply in terms 
of a magnetic scalar potential ip. Because B is also di- 
vergence free, ^ must obey the Laplace equation 



V> = 0. 



(1) 



Since we wish to describe straight beam-line elements, 
it is convenient to work initially in circular cylindrical 
coordinates p, 0, and z with 

x — pcos4), y = psin(f). (2) 

We note for future use that ^ can be written in the form 

x + iy = pe"*' (3) 

from which it follows that, for non-negative integers I and 
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Note that 



CN (z) = d-CtUz)/dz^' 



(9) 



We will require an expansion of ^ in the transverse 
variables x and y. Define the polynomials 



Fl'"" = (x^ + y'^yim{x + iyY 
F^^" = (x^ + y^Ueix + iyY 



(10a) 
(10b) 



for integer I > and m > 0. In Cartesian coordinates we 
then have: 



ip{x,y,z) ^^Yl 



(-l)'m! 

1=0 m=0 ^ ' 

^&^\[z)F'r{x,y) 



'c^\{z)F\^^{x,y) 

(11) 



Note that each polynomial F]f^ is homogeneous of degree 
21 + m. 
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B. Vector Potential 

To determine an associated vector potential A, we 
must find a solution to the coupled system of equations 
V X ^ = VV', where ip is given by the series (111. We 



also must select some particular gauge. And even after a 
particular type of gauge has been chosen, say a Coulomb 
gauge, there is still some remaining freedom [J]. One 
convenient Coulomb gauge choice is given by the rules 



1^0 m^O 

oo oo 
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(12a) 
(12b) 
(12c) 



The coefficients Cm s{z) and Cm dz) describe normal and 
skew components, respectively. For example, in the body 
of a long normal dipole (m = 1), we expect c[°j,(z) 
will be nearly constant (independent of z) and there- 
fore c\^l(z) ~ for n > 0. Correspondingly, in the 



body of a long normal dipole, use of ( 12 1 gives the re- 
sults Aj- ~ Ay c:i and 



A, 



Similarly, in the body of a long normal quadrupole 



2), use of ( 12 1 gives the results A^ — Ay and 



A, 



(13) 



(14) 



For any set of analytic functions Cm,a{z) employed in 



(11), the vector potential defined by (12) satisfies the 
conditions: 1) V x ^ = V^J = B, 2) V x V x A = 0, and 
3) V-j4 = 0. Note that both Maxwell's equations and the 
Coulomb gauge condition are satisfied by construction. 
In the following section, we will see how the coefficient 
bmctions c\^^a{z) can be numerically determined. 



III. 



SURFACE METHODS 



There are cases in which Taylor expansions of the form 



( 11 ) and ( 12 ) can be found analytically. In general, how- 



ever, we have available only measured or numerical three- 
dimensional magnetic field data on a discrete mesh of 
points distributed throughout the region of interest. The 
required high derivatives of V' or A cannot be reliably 
computed directly from this data by numerical differ- 
entiation due to numerical noise, whose effect becomes 
worse with the order of derivative desired. The effect of 
numerical noise, and its amplification by numerical dif- 
ferentiation, can be overcome by fitting on a bounding 



surface far from the axis and interpolating inward us- 
ing the Maxwell equations. Surface fitting methods have 
several advantages, including: 

1. Only functions with known (othonornial) complete- 
ness properties and known (optimal) convergence 
properties are employed. 

2. The Maxwell equations are exactly satisfied. 

3. The results are manifestly analytic in all variables. 

4. The error is globally controlled. Fields that are 
harmonic (fields that satisfy the Laplace equation) 
take their extrema on boundaries. Both the exact 
and computed fields are harmonic. Therefore their 
difference, the error field, is also harmonic, and 
must take its extrema on the boundary. But this 
is precisely where a controlled fit is made. Thus, 
the error on the boundary is controlled, and the 
interior error must be even smaller. 

5. Because harmonic fields take their extrema on 
boundaries, interior values inferred from surface 
data are relatively insensitive to errors/noise in the 
surface data. Put another way, the inverse Lapla- 
cian (Laplace Green function), which relates inte- 
rior data to surface data, is smoothing. It is this 
smoothing that we seek to exploit. In general, the 
sensitivity to noise in the data decreases rapidly (as 
some high inverse power of distance) with increas- 
ing distance from the surface, and this property 
improves the accuracy of the high-order interior 
derivatives needed to compute high-order transfer 
maps. 

Let us briefly compare this approach to that of on-axis 
or midplane fitting. In the case of on-axis fitting, it is 
common to use various analytic model profiles, such as 
Enge functions, and then differentiate them repeatedly 
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to achieve objective 2 by continuing outward using the 
Maxweh equations. However, these functions do not have 
completeness properties, item 1. And there is no smooth- 
ing, item 5, to overcome the amphfication of noise due to 
numerical differentiation. 

In the case of midplane fitting, one approach would 
be to attempt to employ expressions that relate the on- 
axis gradients to midplane data and its derivatives. For 
example, in the case of midplane symmetry, there are the 
relations 



Ct\c{^) = 0, 

cfl{z)^By{x^Q,V = G,z), 
.[0]. ^ 



2 dx 
6 dx^ 



(0,0,; 



(0,0, z) 



1 d^By 
24 5^2 



(15a) 
(15b) 

(15c) 
,etc. (15d) 



{0,0, z) 



By repeatedly differentiating these relations with respect 
to z, one can obtain the Cmjs(-z) for n > 0. In gen- 
eral, the determination of Cm]s{z) requires the computa- 
tion of m -f n — 1 derivatives. Although this approach 
achieves objective 2, since all relevant quantities are sub- 
sequently computed in terms of on-axis gradients, in the 
case where data is available only at grid points it presup- 
poses the ability to compute very high-order derivatives 
by high-order numerical differentiation. This is generally 
impossible, due to the high noise sensitivity associated 
with high-order numerical differentiation, because there 
is no intrinsic smoothing, item 5. 

Another approach is to use an analytic functional form, 
with free parameters, that is known to satisfy the the 3- 
dimensional Laplace equation for all parameter values. 
These parameters can be adjusted so that the field de- 
rived from this representation well approximates the field 
at various grid points. (These grid points could be in 
the midplane, but could be out of the midplane as well.) 
This representation can then be repeatedly differenti- 
ated to provide the required field derivatives. However, 
commonly this fitting procedure has no known complete- 
ness/convergence property, item 1. In some cases Fourier 
series representations with known completeness proper- 
ties are used. But with Fourier series representations, 
an artificial periodicity is imposed in the transverse hor- 
izontal direction. As a result, the Fourier coefficients for 
the field expansions, call them a„, can fall off at best as 
(l/n2)[4]. Correspondingly, the Fourier coeficients in the 
associated expansion for ip fall off at best as (1/n^). As 
a result, repeated differentiation produces nonconvergent 
series, and there is no analyticity, item 3. Whatever rep- 
resentation is used, there is again no intrinsic smoothing 
to overcome the amplification of noise due to numerical 
differentiation. 



A. Use of Field Data on Surface of Circular 
Cylinder 

All three-dimensional electromagnetic codes calcu- 
late all three components of the field on some three- 
dimensional grid. Also, such data is in principle avail- 
able from actual field measurements. In this subsection 
we will describe how to compute the on-axis gradients 
from such field data using the surface of a circular cylin- 
der 13. Once these gradients are known, we may use ( 11 ) 



and ( 12 1 to compute the associated scalar and vector po- 



tentials. 

Consider a circular cylinder of radius i?, centered on 
the 2-axis, fitting within the bore of the beam-line el- 
ement in question, and extending beyond the fringe- 
field regions at the ends of the beam-line element. The 
beam-line element could be any straight element such as 
a solenoid, quadrupole, sextupole, octupole, etc., or it 
could be wiggler with little or no net bending. See Fig. 
[T] which illustrates the case of a wiggler. 




FIG. 1: A circular cylinder of radius R, centered on the z- 
axis, fitting within the bore of a beam-line element, in this 
case a wiggler, and extending beyond the fringe-field regions 
at the ends of the beam-line element. 

Suppose the magnetic field B{x,y,z) — B{p,(f>,z) is 
given on a grid, and this data is interpolated onto the 
surface of the cylinder using values at the grid points 
near the surface. Next, from the values on the surface, 
compute Bp(R, 0, z), the component of B{p, (j), z) normal 
to the surface. The major remaining task is to compute 
the on-axis gradients from a knowledge of Bp{R, (p^ z). 
See Fig. 2. At this point we note that the functions 
exp(ifcz) sin(m0) and exp{ikz) cos(m0) form a complete 
set over the surface of the circular cylinder. 

Let Bp{R, (p, k) be the Fourier transform of Bp{R, (p, z) 
given by the integral 



Bp{R,4>,k) = 



1 

2^ 



^ — ik: 



Bp{R,cp,z). 



(16) 



Because B decays rapidly in the fringe field region, 
Bp{R, 0, z) is absolutely integrable along the z-axis, and 
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Compute Vector 
or Scalar Potential 
Expansion from 
On- Axis Gradients 



Compute On- Axis 
Gradients Using 
Surface Data 
and Inverse 
Laplacian Kernels 



Compute Design 
Trajectory z"^ 
and Transfer 
Map 



FIG. 2: Calculation of realistic design orbit z and its associated realistic transfer map M based on data provided on a 
3-dimensional grid for a real beam-line element. Only a few points on the 3-dimensional grid are shown. In this illustration, 
data from the 3-dimensional grid is interpolated onto the surface of a cylinder with circular cross section, and this surface data 
is then processed to compute the design orbit (trajectory) and the associated transfer map about the design orbit. The use of 
other surfaces is also possible, and may offer various advantages. 
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therefore its Fourier transform is well defined. Next de- 
fine the functions bm,s and bm.c by 



,(i?,fc) 



1 
1 

2^ 



2Tr 



} sm mTT 



) cos TOTT 



dz < 



(17a) 



dz e~^''^Bp{R,(j,,z), 

) 

(17b) 



for m > 1 and 
We know that 



1 

4^ 



27r 



d<t) 



dz e-'^'^^BpiR, (t>,z) 



Bp{R,^,z) = 



dii){p,(j),z) 



dp 



(17c) 



(18) 



p=R 



from which it follows, using the representation ([s]), that 
Bp{R, (t),z)=J2 kl^i{kp)e'''' [G™,,(fc) sinm0 



+Gm,cik) COS mq 



(19) 



Now substitute (191 into the right sides of (17) and per- 



form the indicated integrations to get the results 

bmAR.k) = GrnAk)klUkR), (20) 
from which it follows that 

{R,k) 



Gm,a(k) 



kl'AkR) 



(21) 



This relation for Gm,a{k) can be employed in ^ to give 
the result 



2™to! 



dk 



I'mikR) 



(22) 



We have found an expression for the on-axis gradients 
in terms of field data (normal component) on the surface 
of the cylinder. Equation ( 22 ) may be viewed as the con- 



volution of Fourier surface data b„i^a{R,k) with the in- 
verse Laplacian kernel k"^"^~^ /I'AkR)- Moreover, this 
kernel has a very desirable property. The Bessel functions 
I'AkR) have the asymptotic behavior 



\I'AkR)\ - exp(|fc|i?)/v/27r|fc|i?as |fc| 



— oo. 



(23) 



Since I'^kR) appears in the denominator of ( p2[ ), we see 
that the integrand is exponentially damped for large 
Now suppose there is uncorrelated point-to-point noise in 
the surface data. Such noise will result in anomalously 
large |fc| contributions to the bm,a{R, k). But, because 
of the exponential damping arising from I'^kR) in the 
denominator, the effect of this noise is effectively filtered 



out. Moreover, this filtering action is improved by mak- 
ing R as large as possible. This filtering, or smoothing, 
feature will be discussed in more detail in Section V. 

We close this subsection with the remark that if one 
wishes to extract the Cq^I{z) (monopole) on-axis gradi- 
ents from field data, as is required for example in the case 
of a solenoid, it may be preferable to use the longitudi- 
nal component Bz{R, (j), z) on the surface of the cylinder 
rather than the normal component Bp{R, 4>,z) 



B. Use of Field Data on Surface of Elliptical 
Cylinder 

1. Background 

In the previous subsection we employed a cylinder with 
circular cross section, and observed mathematically that 
it is desirable for error insensitivity to use a cylinder with 
a large radius R. Physically, this is because we want the 
field data points employed to be as far from the axis 
as possible since the effect of inhomogeneities (noise) in 
the data decays with distance from the inhomogeneity. 
Evidently the use of a large circular cylinder is optimal 
for beam-line elements with a circular bore. However, 
for dipoles or wigglers with small gaps and wide pole 
faces, use of a cylinder with elliptical cross section should 
give improved error insensitivity. See Fig. [3] In this 
subsection we will set up the machinery required for the 
use of elliptical cylinders, and apply it to the calculation 
of on-axis gradients based on field data [3J|3j. 




FIG. 3: An elliptical cylinder, centered on the z-axis, fitting 
within the bore of a wiggler, and extending beyond the fringe- 
field regions at the ends of the wiggler. 



2. Elliptic Coordinates 

Elliptic coordinates in the x,y plane are described by 
the relations [Sj 



X = f cosh(u) cos(z;), 
y = /sinh(M) sin(t!). 



(24a) 
(24b) 



Contours of constant u, with u € [0,oo], are nested el- 
lipses with common foci located at {x\ y) = (±/; 0). Con- 
tours of constant with v £ [0, 27r], are hyperbolae. 
Together these contours form an orthogonal coordinate 
system. See Fig. [4j Data is to be interpolated onto the 
ellipse whose cross section is that of the elliptical cylinder 
of Fig. |3l See Fig. (5] 




FIG. 4: Elliptical coordinates showing contours of constant 
u and constant v [7]. The foci are at (±a, 0), and in our case 
a = f. 




FIG. 5: A square or rectangular grid in the x,y plane for a 
fixed 2 value on the 3-dimensional grid. Values at data points 
near the ellipse are to be interpolated onto the ellipse. 

For our work we will need the unit vector e^, the unit 
vector (outwardly) normal to the surface of the elliptical 
cylinder. Write 



yey 



f cosh(M) cos(i')e2: + / sinh(u) sm{v)ey + ze^ 



Then, by definition, we have the result 

e„ = {dr/du)/\\{dr/du)\\ 

sinh(u) cos(w)e^ -I- cosh(M) sin(w)ej, 
[cosh^(M) — cos^(u)]i/2 



(25) 



(26) 



It is also convenient to employ the complex variables 
C = x + iy, (27a) 

and 

w — u + iv. (27b) 



In these variables, the relations (24 1 can be written in 
the more compact form 



C — /cosh(w). 



(28) 



Form differentials of both sides of ( |28[ ). Doing so gives 
the result 

dx + idy — f sm]i{w) {du + idv) (29) 
and the complex conjugate result 

dx — idy = f sinh('i())(du — idv). (30) 



Now form the product of ( 29 ) and ( |30| ) to get the trans- 
verse line-element relation 



dsi 



dx^ + dy^ — sinh(u -f iv) sinh(u — iv){dv? + dv"^) 



= f[cos\i^{u) - cos2(u)](du2 + dv'^). 
From this relation we infer the results 



(31) 



= (l//)[cosli2(u) - cos^{v)]-^'^{di^/du), (32a) 

= (l/,f )[cosh2(w) - cos2(«)]-i[(a„)2 + (a,)2]V' 
+ {d.)^i^- (32b) 

3. Mathieu Equations 

Let US seek to construct harmonic functions of the form 

iP P{u)Q{v)exp{ikz) (33) 

where the functions P and Q are yet to be determined. 
Employing the Ansatz (33 1 in Laplace's equation and use 



of ( 32b ) yields the requirement 

[(du)' + {d,)'mu)Q{v)] = 

fcV^[cosh^(u) - cos^{v)]P{u)Q{v). (34) 

We also observe that there is the trigonometric identity 

cosh^ (u) " cos^ (v) = (1/2) [cosh(2u) - cos(2w)] (35) 

so that the requirement ( |34[ ) can be rewritten in the form 

P„)2 + (9„)2][P(u)Q(t;)] = (36) 
{Pf/4)[2 cosh(2u) - 2cos(2w)]P(u)Q(i;). 
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Upon dividing both sides by PQ, (36) becomes 

{l/P){d^fP+{l/Q){d,YQ^ 

(fc2/V4)[2cosh(2u) - 2cos(2i;)], (37) 

from which it follows that 

(1/P)(a„)2p_ (fc2j2/4)[2cosh(2?/)] = (38) 
- (1/Q)(9„)2g- (fc2/V4)[2cos(2z;)]. 

Therefore, there is a common separation constant a such 
that 



and 



(l/P)(a„)^P- (fcV74)[2cosh(2u)] = a (39a) 



(l/Q)(a„)2Q- (fc2/2/4)[2cos(2w)] =a. (39b) 



Correspondingly, P and Q must satisfy the ordinary lin- 
ear differential equations 

cPP/du^ -[a-2q cosh(2u)]P = 0, (40a) 



where 



cfQ/dv^ + [a-2q cos{2v)]Q = 0, (40b) 



k'^f/A. (40c) 



Equation (40b) for Q is called the Mathieu equation, 
and Equation (40a) for P is called the modified Mathieu 



equation. For our purposes, we will need solutions Q{v) 
of (40b) that are periodic with period 27r. Such solutions 



exist only for certain characteristic values of the sepa- 
ration constant a. These values are denoted a„(q) for 
n = 0, 1, 2, 3, • • • and 6„(g) for n = 1, 2, 3, • • • . The solu- 
tions associated with the separation constants a = an{q) 
are denoted ce„(w, q). They are even functions of v and, 
in the small q limit, are proportional to the functions 
cos{nv). The solutions associated with the separation 
constants a — 6,1 (?) are denoted se„ (11,(7). They are odd 
functions of v and, in the small q limit, are proportional 
to the functions sin(nu). The functions ce„(t;,g) and 
se„(w, q) form a complete orthogonal set over the interval 
V G [0, 27r] and are normalized so that 



dv ce„(u, g) ce„(w,q) = 7r(5„j„, (41a) 



2tt 



dv se„(u, q) se„(w, q) = irSmn, (41b) 



2lT 



dv ce„(u, q) se„(i;, q) = 0. (41c) 



With regard to the solutions of the modified Mathieu 



equation, note that ( 40b ) is transformed into ( 40a I under 
V — ^ iu. As a result, corresponding (real-valued) solu- 
tions to (40a) are defined by Ce„(M, g) = ce„(jM, g) and 
Se„(M,g) = —isem{iu,q). We refer the reader to O |4] 
and [5] for a detailed treatment of the Mathieu functions 
and their properties. 



4- Elliptic Cylinder Harmonic Expansion and On- Axis 
Gradients 

The stage is now set to describe the expansion of any 
harmonic function tjj in terms of Mathieu functions. The 
general harmonic function that is analytic in x and y near 
the origin can be written in the coordinates ( 24 ) in the 
form 



ip{u, V, z) ~ 



n=0 

0° poo 

E 

n=l 



dk c„(fc)e*''''Ce„(u, q) ce„(w, q) 
dk s„(fc)e*''^Se„(u, 5) sen{v,q) 



(42) 

where the functions Cn{k) and Sn(fc) are arbitrary. We 
will call (42) an elliptic cylinder harmonic expansion. 



To exploit this expansion, suppose the magnetic field 
B{x,y,z) is interpolated onto the surface u = U oi 
an elliptic cylinder using values at the grid points near 
the surface. See Fig. [5j Let us employ the notation 
B{x, y, z) = B{u, V, z) so that the magnetic field on the 
surface can be written as B{U,v, z). Next, from the 
values on the surface, compute Bu{U,v, z), the compo- 
nent of B(x,y,z) normal to the surface. Our aim will 
be to determine the on-axis gradients from a knowledge 
of Bu(U,v, z). At this point we note that the functions 
exp(ifcz)se„(f , g) and exp(ifcz)ce„(w, g) form a complete 
set over the surface of the elliptical cylinder. 

Let us begin by solving (32a I for {d^/du). We find, 
using ([26]), the result, 

{dil)/du) = f[cos\?{u)-cof?{v)f'^Bu (43) 
— /(sinh u cos i')i?a; + /(cosh M sin w)i?j,. 



We see that the right side of ( 43 ) is a well-behaved func- 
tion -F(m, f , z) whose values are known for u — U, 

F{U,v,z) — /(sinh U cos v)Bx {U, v, z) 

+ f{coshUsinv)By{U,v,z). (44) 



Moreover, using the representation (^42h in (43) and 



( 44 ) , we may also write 



FiU,v,z) 



n=0 



dk Cn{k)e'''''Ce'„{U, q) ce„(w, q) 



00 pQQ 

dk s„(/c)e*'=^SeUC/, q) se„(z;, q). 

1 -/ —00 



(45) 



Next multiply both sides of (45) by exp{—ik'z) and inte- 
grate over z. So doing gives the result 

— / dz e-'^'F{U, «, z) = V c„(fc)Ce;(t/, q) ce„(«, q) 

00 

+ ^s„(fc)Se; ([/,(/) 

n=l 

(46) 
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Now, employ the orthogonality properties of the Mathieu 
functions to obtain the relations 



Crik)Ce'^iU,q) = 



27r2 



27r 



dvcer{v,q) / dz e ''^'^ F(U,v, z), 



(47a) 



s,(fc)Se;(f/,(7) = 



27r2 



27r 



V, z ) 



(47b) 



In view of (471, define the function F{v^ k) by the rule 



F{v,k) = 



1 
2n 



dz e-'^'F{U,v,z), (48) 



and define functions F^{k) and F^{k) by the rules 

r2n 



1 



27r2 



dv cer{v, q)F{v, k) 



(49a) 



dv ser{v, q)F{v, k) 



27r2 



27r 



dvser{v,q) / dz e ''''^ F{U,v, z). 



(49b) 



Note the similarity to ( |17[ ), where cos(r0) and sin(r0) 
are replaced by cer{v,q) and ser{v,q). We will call the 

functions F" (k) Mathieu coefficient functions in analogy 
to the Fourier coefficients that arise in Fourier analysis. 



With these definitions, the relations (47 1 can be rewrit- 
ten in the form 



c (k) - . (k) 



(50) 



Finally, employ ( 50 ) in ( 42 ) . So doing gives the result 



J2 / dk e'''''[F^{k)/Ce'^{U,q)]Cer{u,q) cer{v,q) 

+ J2 dke''''[F^{k)/Se'^{U,q)]Ser{u,q)ser{v,q). 

(51) 

We have obtained an elliptical cylinder harmonic expan- 
sion for -ijj in terms of surface field data. 

Of course, what we really want are the on-axis gradi- 
ents. Again, once these gradients are known, we may use 



(111 and (12 1 to compute the associated scalar and vec- 



tor potentials. The gradients can be found by employing 
two remarkable connections (identities) between elliptic 
and circular cylinder functions of the form [3] 



oo 



Cer{u,q) cer{v,q) — > a^^{k)Im{kp) cos.(m(j)), (52a) 



?n=0 



Ser(u,q) ser{v,q) = ^ I3l^{k)l.m{kp) sin{m(t)). (52b) 



For further reference, we will call the quantities a^(A;) 
and Pm{k) Mathieu-Bessel connection coefficients [3l|4]. 
Using these results, (51) can be rewritten in the form 



^p{x,y,z) = 

J2 dk e'^'Irnikp) cos(m0) ^ a;,(fc)[F;(fc)/Ce;(C/, q)] 



m=0 • 

oo „oo 



dv cer{v,q) I dz e F{U,v,z), rn=i 



r=Q 
oo 



J2 I dk e'^'Irnikp) sm{m(b) Pl,{k)[K i^) I KiU , q)]- 

" r=l 

(53) 



Upon comparing (53) with ([5]), we conclude that there 
are the relations 



G™,e(fc) ^Y.<^^^)^^r{k)IGe'AU,q% (54a) 



and 



GmAk) = ^/3;,(fc)[#;(fc)/Se;(C/,g)]. (54b) 



Finally, in view of (H and ([54]), we have the desired re- 
sults 

CtU^) = 

■n fcc oo 

— / *e^'=^fc"+"^a:„(fc)[F;(fc)/Ce;(f/,g)], 

^ "i- ^-oo 

(55a) 

■n roo oo 

— / dfce''=^fc"+"^/?:„(A:)[F,^(fc)/Se;(L/,g)]. 

(55b) 

We have found expressions for the on-axis gradients in 
terms of field data (normal component) on the surface of 
an elliptic cylinder. 
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C. Use of Field Data on Surface of Rectangular 
Cylinder 

A similar procedure has been developed for computing 
the on-axis gradients in terms of field data provided on 
the surface of a rectangular cylinder. In this case, each 
on-axis gradient Cm]a may be written as the sum of four 
contributions, 

Ct]M)^ J2 (56) 

0=T,B,LM 

Each contribution is determined by the integration of the 
normal component of the field against an appropriate ker- 
nel over one of the four faces (Top, Bottom, Left, Right) 
of the rectangular cylinder. The resulting expressions are 
quite lengthy, and we therefore refer the reader to O |4] 
for further details. The remainder of this paper will focus 
on the circular and elliptical cylinder cases. 



y 



+g ' 


X 


/ — 


/I 




u 


-g ' 





FIG. 6: A monopole doublet consisting of two magnetic 
monopoles of equal and opposite sign placed on the y axis 
and centered on the origin. Also shown, for future reference, 
is a cylinder with circular cross section placed in the interior 
field. 



IV. NUMERICAL BENCHMARKS 
A. Monopole Doublet 

In this section, we develop an exactly-soluble but nu- 
merically challenging model field to be used to numeri- 
cally benchmark the procedures developed in Section III. 
Suppose two magnetic monopoles having strengths ±g 
are placed at the (cc, y, z) locations 



+g 



r+ = (0,a,0), 
r- = (0,-a,0). 



(57) 



See Fig. [6j which also shows a circular cylinder with 
radius R (the surface p — R). These monopoles generate 
a scalar potential "0(2^1 ?/? described by the relation 

il){x,y,z) = 

-g[x^ + {y - «)' + ^^r"^ + 9[^^ + {V + «)' + ^^"^ 



(58) 



Correspondingly, they produce a magnetic field B — "S/il^ 
having the components 

= gx[x'^ + {y-af + z'^]-^/'^-gx[x^ + {y+af+z^Y^''^, 

(59a) 

By=g{y-a){[x^ + {y~af + zX^'^ 

~g{y + a)[x^ + {y + af + z^r^'^}, 

(59b) 



= gz[x^ + {y-af + z^]-''/^-gz[x^ + {y+af + z^]-^'\ 

(59c) 

This field is sketched in Fig. [7| To provide fur- 
ther insight. Fig. [8] shows the on-axis field component 
By{x = 0,y = 0,2), and Figs. [9] and 10 show the 




FIG. 7: The interior field of a monopole doublet. Also shown 
is an eUipse which will be used in Section [IV C[ 



off-axis field components Ba:{p — 1/2,0 — 7r/4, z) and 
B,(p= 1/2,0 = 7r/4,z). 

Due to the symmetries of the field, the functions Cm]c 
vanish for all m. Furthermore, Cm]s = for m even. It 
can be shown that the nonvanishing on-axis gradients are 
given for m odd by [3]: 

rM (~\-( 9 (2™)! n2m+l(-. 

(60a) 

where 



/3(2) 



Vz^ + 



(60b) 



For fixed z, the domain of convergence for the poly- 
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By(T) 




FIG. 8: The on-axis field component By{x = 0,y = 0, z) for 
the monopole doublet in the case that a = 2.5 cm and g = 1 
Tesla-(cm)^. The coordinate z is given in centimeters. 



Bx (T) 




z (cm) 



FIG. 9: The field component on the line p = 1/2, <^ = 7r/4, 
z G [~oo, oo] for the monopole doublet in the case that a — 2.5 
cm and g = 1 Tesla-(cm)^. The coordinate z is given in 
centimeters. 



Bz (T) 
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FIG. 10; The field component on the line p — 1/2, (f> = 
7r/4, z £ [—00,00] for the monopole doublet in the case that 

\2 



a = 2.5 cm and g ~ 1 Tesla-(cm) 
in centimeters. 



The coordinate z is given 



nomial series ( H][T2 I representing the field in terms of 



V z (cm) / 



the functions ( 60 1 is given by the condition i/a;^ + j/^ < 



+a^. In particular, the domain of convergence is 
a region of circular cross-section whose radius increases 
as we move longitudinally away from the location of the 
monopoles at z = 0. 

Suppose we wish to calculate a transfer map through 
jth Qj-^Qj-. Then, as shown in the Appendix, to do so 
requires knowledge of the Cm]a{z) with {m + n) < 7 

when m = or m is odd, and knowledge of the Cm]a{z) 
with (m + n) < 8 when m is even. 

Graphs of a selected few of these functions, for the 
monopole doublet in the case th at a = 2.5 c m a nd 5 = 1 
Tesla-(cm)^, are shown in Figs. 
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through 
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In these 

plots z has units of centimeters. Evidently the Cm.s be- 
come ever more highly peaked with increasing m. For- 
tunately, when working through some fixed degree, we 
need fewer derivatives with increasing m. Note that we 
expect that the function Cm]siz) should have n zeroes. 
This is indeed the case. 



B. Circular Cylinder Results 



The procedure discussed in Section III A has been 
benchmarked using the field of a monopole doublet in 
the case that a = 2.5 cm and g = 1 Tesla-(cm)^. We 
set up a regular grid in x, y, z space, where we let 
each variable range over the intervals x G [—4.4, 4.4] 
with spacing = 0.1, y € [—2.4,2.4] with hy = 0.1, 
and z £ [—300,300] with = 0.125 (in units of cm). 
The known values of the three components of the field 
are computed using (59) at each grid point. Consider 
a cylinder of radius R = 2 cm and length 600 cm. We 
use bicubic interpolation to interpolate B at these grid 
points onto 49 selected angular points on the cylinder, 
for each of the 4801 selected values of z. The angular in- 



tegration in (17) is performed using a Riemann sum with 
iV = 49. (This is necessary to ensure sufficient conver- 
gence of the angular integrals to within lO"**.) We evalu- 
ate the Fourier transform at 401 values of k in the range 
[—Kc, Kc] with Kc = 20, using a spline-based Fourier 
transform algorithm [31 0] . We use these same points in 
k space to evaluate the inverse Fourier transform, provid- 
ing a set of numerically determined functions Cm]a{z). 



A comparison between the exact on-axis gradients ( 60 ) 



and th ose obtained from grid data is provided in Figs. 

Evidently 
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for the functions cfl, C|"|,, and Cy 
the agreement is excellent. Fig. [14| illustrates the differ- 
ence between the on-axis gradient c[°j, as obtained from 
grid data and the exact on-axis gradient obtained from 
(60) with m = 1. The maximum error attained rela- 
tive to peak is 1.7 x 10~^. Further detailed study shows 
that numerical results agree with exact results to within 
relative errors less than a few parts in 10^ for all the 
relevant C'm]a{z). By using exact data on the cylinder 
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rather than interpolating off the grid onto the cyhnder, 
we have also verified that the error due to interpolation 
onto the cylinder is comparable to that produced by nu- 
merical integration [3^. Finally, all these small errors can 
be further reduced with the aid of a finer grid H] . 
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FIG. 13: Exact and numerical results for C^j\{z). Exact 
results are shown as a solid line, and numerical results are 
shown as dots. 



FIG. 11; Exact and numerical results for C^^\{z). Exact 
results are shown as a solid line, and numerical results are 
shown as dots. 





C 4.7 X 10'^ 
O 4.6 X 10'^ 
LU 4.5X10"^ 



FIG. 14: Difference between exact and numerical results for 



FIG. 12; Exact and numerical results for c{''!,(z). Exact 
results are shown as a sofid fine, and numericaf resufts are 
shown as dots. 



C. Elliptical Cylinder Results 



The procedure discussed in Section |III B| has been 
benchmarked using grid values identitical to those de- 
scribed in the previous section. Consider an elliptical 
cylinder of semiminor axis of ^max — 2 cm and semima- 
jor axis Xmax — 4 cm. In this case, we evaluate the an- 



gular integrals (49) using a Riemann sum with N — 120. 
(This is necessary to ensure sufficient convergence of the 
angular integrals to within 10~*.) Doing so requires in- 
terpolation off the grid onto the elliptical cylinder at 120 



angular points for each of the 4801 selected values of z. 
The sums in ( 55 ) are truncated beyond r = rmax , where 
fmax varies from 11 to 29 as necessary to achieve a toler- 
ance of 1 part in 10'*. We evaluate the Fourier transform 
at 401 values of k in the range [—KcKc] with Kc — 20, 
using a spline-based Fourier transform algorithm. We 
use these same points in k space to evaluate the inverse 
Fourier transform, providing a set of numerically deter- 
mined functions C^^a{z)- 

Results for the functions Cm|Q(z) are similar to those 
found in the circular cylinder case [31 S] , and have com- 
parable accuracy. In the following section, however, we 
illustrate that functions obtained using an elliptical cylin- 
der are significantly more robust against numerical noise 
in the original grid values. 
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V. SMOOTHING 

In this section, we investigate the smoothing of numer- 
ical noise that results from the use of the surface fitting 
techniques described in Section Consider computing 
the on-axis gradients C^m,a from a grid of numerical field 
values B using a circular or elliptical cylinder, as de- 
scribed in the previous sections. Observe that (22 1 and 



according to (61). Figs. 19p0 illustrate the on-axis gra- 



( 55 1 are linear in each of the values Bj. and By at the grid 



points. The inclusion of additive numerical noise AB at 
each grid point therefore results in on-axis gradients of 
the form Cmja -I- ACmJa, where the contribution ACmL 
is determined by the values AS according to the proce- 
dures of Section [Hi) Note that only the field values AB^ 
and ABy at grid points near the surface of the circular 

cylinder affect the functions ACm|a. 

To examine the effect of additive noise, we gener- 
ate a random noise field AB whose components are 
proportional, at the 1% level, to the strength of the 
monopole-doublet on-axis vertical field. Consider the 
grid [-4.4, 4.4] x [-2.4, 2.4] x [-300, 300] cm used in Sec- 



tion IV for fitting the field of the monopole doublet, with 
mesh points indexed by j = 1, • • • A^. Let i?y(0, 0, z) de- 
note the value of the monopole-doublet on-axis vertical 
field at longitudinal location z, as determined from (59). 
At each point [xj, j/j, Zj) we set 

ABx{xj,yj,Zj) ^ eBy{0,0,Zj)Sx{j), (61a) 
ABy{xj,yj,Zj) = eBy{0,0,Zj)Syij). (61b) 

Here the Sx{j) and dy(j) are uniformly distributed ran- 
dom variables taking values in the interval [—1, 1], and 
e — 0.01. After interpolating these values onto the sur- 
face of a circular cylinder with R = 2 cm and z e 
[—300, 300] cm, we use the procedure described in Section 



III A to compute the on-axis gradients (22). 

In Fig 
tity bm,s 



|15[ we have displ aved the computed quan- 
s{R,k) appearing in (17) for the case m = 1. 
It is a function of the spatial frequency fc, having ran- 
dom variations of approximately uniform variance over 

-1 



16 displays the kernel 



the interval [-20,20] cm~\ Fig. 

fc'"-V/;„(fci?) multiplying 6,„,,(i?,fc) in ^ for the case 
m — 1 and R = 2 cm. Note the rapid decay of this func- 



tion for large |fc|. Finally, Fig. 17 displays the product 
of these functions, illustrating the dramatic suppression 
of high-fc contributions to the Fourier integral appearing 
in ((22l). 



A similar phenomenon occurs when fitting is performed 
using an elliptical cylinder. In this case, a sequence of 



kernels contributes to (55) for each fixed m. In Fig. 18 



we have displayed the kernels contributing to the case 
m — \, with Xmax = 4 cm and ^max = 2 cm. Kernels 
take their maxima at A: = 0, and these maxima decrease 
monotonically with increasing index r. All kernels de- 
crease rapidly with increasing |fc| [31 E]. 

To study the effect of noise on the on-axis gradients, 
twelve distinct random fields were generated on a mesh 



dients C^), and as computed using these field values. 
The solid line in Fig. [19] illustrates the rms value of the 



on-axis gradient C{ 



[6] 



as computed using a circular cylin- 
The dashed 



der of radius R = 2 cm according to (|22 ) 
line in Fig, 



gradient cf'],, as computed using an elliptical cylinder of 
semiminor axis 2 cm and semimajor axis 4 cm according 
to ([55l). 



19 illustrates the rms value of the on-axis 



In Fig. 



similar results are shown for the 

on-axis gradient C}"^ 

The attentive reader might wonder why we have dis- 
played the Cm,c for noise while the the monopole doublet 

In] 

field is governed by the gradients Cm 's ■ The reason is that 

In] 

the Cm's are produced by fields that are predominantly in 
the vertical (y) direction, and for such fields there is rel- 
atively less difference between circular and elliptical sur- 
face fitting because the semi-minor axis of the elliptical 
cylinder is the same as the radius of the circular cylinder. 
Note also that in both cases only the component of the 
field normal to the surface is used. Horizontal fields drive 
primarily the Cm]c- However, in the horizontal direction, 
the semi-major axis of the elliptical cylinder is substan- 
tially larger than the radius of the circular cylinder. We 
therefore expect the advantage of using elliptical cylin- 
ders compared to circular cylinders will be most apparent 

\n] 

in the Cm'c- Indeed, this is what Figs. 19 and 20 show. 
The effects of errors in the surface data are suppressed 
more when the bounding surface is farther away from the 
field observation point. As a result of this suppression, 
it is advantageous to use a fitting surface that is as far 
away as possible. 

Suppose we compare the Cm]c due to noise and shown 

In] 

in Figs. 19 and 20 with the corresponding Cm's in 
Figs. 12 and 13. This is reasonable because a horizontal 
monopole doublet would produce Cm]c analogous to the 
Cm]s shown in Figs. 12 and 13. Remarkably, we find 
that, as a result of smoothing, the effect of 1% noise in 
the field data produces, on average, only on the order of 
0.01% changes in the on-axis gradients. 



VI. APPLICATIONS 

A. ILC Damping Ring Wiggler Fields 

A less stringent test of the accuracy of surface meth- 
ods is that the magnetic field, as computed from surface 
data using field values on a 3-dimensional mesh, should 
reproduce the magnetic field at the interior mesh points. 
(This is also a test of the quality of the magnetic data 
on the mesh.) We computed such an interior fit, and the 
associated transfer map, for the modified CESR-c design 
of the Cornell wiggler, which has been adopted as the 
design prototype for use in International Linear Collider 
studies [HI HO] ■ Cornell provided data obtained from the 
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k(cm"') 




FIG. 15: The quantity 6i,s(i?, fc) computed from uniform 
random noise (61 1 using a circular cylinder with R = 2 cm. 
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FIG. 16: The kernel k"^~^ / I'„,{kR) as a function of k for the 
case m = 1 and R = 2 cm. 
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FIG. 17: The product k"^-'^bn^,4R,k)/I',„ appearing in (22l 



as computed from uniform random noise ( 61 1 with m = 1 and 
i? = 2 cm. 



FIG. 18: The kernels k"' l3;^{k) /Se'r{U, q) for the case m = 1 
and r = 1, 3, 5, 7, 9, 11, as a function of k, with q and k related 
by (liOcf. 
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FIG. 19: The function C[ I as computed from a mesh contain- 
ing random noise (61 1 at each mesh point. Damping of this 



noise illustrates the effect of smoothing. (Solid line) Result 
obtained using a circular cylinder with R = 2 cm. (Dashed 
line) Result obtained using an elliptical cylinder with semimi- 
nor axis of 2 cm and a semimajor axis of 4 cm. 



3-diniensional finite element modeling code 0PERA-3d 
for the field components B^, By, and B^ on a mesh of 
spacing 0.4 x 0.2 x 0.2 cm in a volume 10.4 x 5.2 x 480 
cm, extending beyond the fringe-field region. The field 
components are provided to a precision of 0.1 G rela- 
tive to a peak field of 16.7 kG. An elliptic cylinder with 
semimajor axis 4.4 cm and semiminor axis 2.4 cm was 
placed in the domain of the data, and the field on the 
elliptic cylinder boundary was constructed using bicubic 
interpolation. See Figs. [3] and [5] 

The interior field was computed using the on-axis gra- 
dients through terms of degree 6 in x,y over the do- 
main of the original data. This solution for the interior 
field was then compared to the original data at each grid 



point. Fig. 21 displays the vertical field By off-axis at 
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FIG. 20: The function C| I as computed from a mesh contain- 
ing random noise (61 1 at each mesh point. Damping of this 



noise illustrates the effect of smoothing. (Solid line) Result 
obtained using a circular cylinder with R = 2 cm. (Dashed 
line) Result obtained using an elliptical cylinder with semimi- 
nor axis of 2 cm and a semimajor axis of 4 cm. 



{x,y) = (0.4,0.2) cm along the length of the wiggler. 
The field data (points) are shown along with computed 
values (solid line) . Note that the fitted field captures the 
fringe-field behavior. The relative error was found to sat- 
isfy the bound (5|Bdata - B/»t|/|-B|peafe < 3.5 x 10-^. We 

observe that this error is comparable to that found for 
the monopole-doublet benchmark. Presumably it is due 
to errors in numerical integration, errors in interpolating 
onto the elliptic cylinder, errors arising from neglecting 
terms beyond degree 6, etc., as well as possible failure 
of the 0PERA-3d data to be Maxwelhan. Fig. [22] illus- 
trates the horizontal roll-off of the vertical field at y = 0.1 
cm, z — 104.2 cm for < a; < 1 cm. Note the discrete 
jumps in the original data, reflecting the number of digits 
retained in the output of the numerical computation. De- 
spite the small variation of By in x, the fit goes through 
the interior data. Finally, Fig. [23] illustrates the longi- 
tudinal field Bz, again at (x, y) = (0.4, 0.2) cm along the 
wiggler. Note that no information about Bz was used to 
generate this field, since only the component of B normal 
to the elliptic cylinder surface was used to generate the 
interior solution. 

The error for By on-axis lies in the range 0.1-0.2 G 
along the length of the wiggler, increasing slightly near 
the end poles. A plot of residuals in the plane y — is 
displayed in Fig. [24] Note that the error is within 0.6 
G over this region of the x-z plane. The error begins 
to increase rapidly at about x = 2 cm; this may be due 
to retaining only terms through degree 6 in the on-axis 
gradients, or perhaps a finite domain of convergence of 
the associated power series for By{x, y, z). For all \x\ < 2 
cm, the peak error is 0.3 G. We remark that this peak 
error amounts to a relative error of less than 2 parts in 
10^, which is remarkably small compared to the error 
for the data of Fig. 21. We expect the error to behave 
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FIG. 21: Fit to the proposed ILC wiggler vertical field versus 
z, where x = 0.4 cm, y = 0.2 cm. The solid line is computed 
using data on the surface of an elliptical cylinder with Xmax = 
4.4 cm, ymax ~ 2.4 cm, using the polynomial series for B 
obtained from (111 or (12 1. Dots represent numerical data 



provided by OPERA-3d. 
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FIG. 22: Fit to the proposed ILC wiggler vertical field versus 
X, where y = 0.1 cm, z — 104.2 cm. The solid line is computed 
using data on the surface of an elliptical cylinder with x,nax = 
4.4 cm, ymax = 2.4 cm, using the polynomial series for B 



obtained from (111 or (12 1. Dots represent numerical data 



provided by OPERA-3d. 



like a harmonic function, and therefore it should grow 
as one approaches the boundary. Conversely, it should 
be the smallest on the center line. The observed error 
follows this pattern. Finally, this phenomenon may also 
be a factor in the observed increase in the error at and 
beyond x = 2 cm. 
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FIG. 23; Fit to the proposed ILC wiggler longitudinal field 
versus z, where x — 0.4 cm, y = 0.2 cm. The solid line is com- 
puted using data on the surface of an elliptical cylinder with 
Xmax = 4.4 cm, ymax = 2.4 Cm, using the polynomial series 
for B obtained from or ( |12[ |. Dots represent numerical 
data provided by 0PERA-3d. 



Difference in By (Gauss) 





FIG. 24: Difference (Gauss) between the vertical field By 
of the proposed ILC wiggler and its fitted value across the 
midplane y = 0. Peak field is 16.7 kG. 



B. ILC Damping Ring Wiggler Design Orbit and 
Associated Transfer Map 

The on-axis gradients computed for the ILC wiggler 
were then used in the code MaryLie/IMPACT to inte- 
grate, simultaneously, i) equations for the design orbit of 
a 5 GeV positron through the wiggler, ii) equations for 
the matrix elements of the linear part of the transfer map 
through the wiggler, and iii) equations for the coefficients 
of the generating polynomials /3,/4,... appearing in the 



-0.2 

X (mm) 



FIG. 25: Design orbit for a 5 GeV positron through the 
proposed ILC wiggler. (Upper) Coordinate a:;(mm) along the 
length of the wiggler z{m). (Lower) Design orbit in the phase 
space defined by coordinates a;(mm) and Px/p^- 



Lie factorization of the transfer map. Each generator of 
the symplectic transfer map M is computed in variables 
representing deviations from the design orbit. See the 
Appendix for a brief discussion of Lie methods. 

The design orbit is displayed in Fig. [25j and Table I 
lists its initial and final conditions. Table II displays the 
matrix i?2 that describes the linear part of the transfer 
map in (62). Phase-space coordinates are arranged in 
the order (x, Pa;, y,Py, t,Pt). The first few Lie generators 
/,„ of the nonlinear part of the transfer map are listed in 
Table III. 



VII. 



CONCLUDING SUMMARY 



Surface methods provide a reliable and numerically ro- 
bust technique for extracting transfer maps from numer- 
ical field data. By benchmarking them against a nu- 
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TABLE I: Design Orbit Specifications. 



E [GeV] 


[GeV/c] 


Bp [Tm] 


5.00000000000000 


5.000510972950664 


16.67990918220719 




Entry {z = m) 


Exit (2 = 4.8 m) 


X [cm] 


0.0 


-0.00449920546655753 




0.0 


4.5494526506306504x10"^" 


V [cm] 


0.0 


0.0 




0.0 


0.0 


ct [cm] 


0.0 


480.00494875303383 




-1.0000000052213336 


-1.0000000052213336 



TABLE IL Linear transfer map R2 for the ILC damping ring wiggier. 

1.000056 4.800233 x 10^ 0.000000 0.000000 0.000000 -4.500058 x 10"^ 

2.235501 X 10"^ 1.000052 0.000000 0.000000 0.000000 -1.001075 x 10"^ 

0.000000 0.000000 9.404373 x 10"^ 4.687866 x 10^ 0.000000 0.000000 

0.000000 0.000000 -2.465383 x lO"* 9.404414 x 10"^ 0.000000 0.000000 

-4.857478 x 10"" -4.499810 x 10"^ 0.000000 0.000000 1.000000 9.897806 x 10"^ 

0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 



TABLE III: First few nonvanishing Lie generators fm for the 
ILC damping ring wiggier. 

Index Monomial Coefficient 
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-1.0738513995490168X 10"^ 


29 


X^Px 


4.630793805976143x10-'' 


33 


X^Pt 


1.1070545309022173x10"^ 


34 


xpl 


-2.2214155912078926X 10"^ 


209 


pt 


-4.949353522256519 x lO^^* 



merically challenging problem whose results are known 
exactly, we have verified that surface methods have all 
the advantages claimed in the beginning of Section III. 
In particular, we demonstrated that errors of a few parts 
in lO"*^ can be achieved for all on-axis gradients Cm]a{z) 
required to compute transfer maps through 7*'^ order, 
and that the results obtained were remarkably insensi- 
tive to noise. Moreover, these small errors can be further 
reduced if desired with the aid of a finer grid. 

Subsequently we applied surface methods to compute 
the interior field for the proposed ILC Damping Ring 
wigglers. Consistent with the accuracy displayed by the 
monopole-doublet benchmark results, excellent fits were 
demonstrated for interior fields. We also illustrated the 
computation of the design orbit and its associated trans- 
fer map based on surface methods. 

In summary, the use of surface methods makes it pos- 
sible, for the first time, to compute for straight beam-line 
elements realistic transfer maps for real magnets includ- 



ing all fringe and high-order niultipole field effects. 

In many cases, however, we are interested in magnetic 
elements with significant sagitta, such as dipoles with 
large bending angles. In these cases, it is not possible 
in general to surround the design orbit with a cylindri- 
cal surface that lies interior to all iron or other magnetic 
sources. Part II of this paper will describe an alterna- 
tive, but more computationally intensive, method suit- 
able for general geometries including that of a bent box 
with straight ends. See Fig. [26] In this case we obtain 
simple, geometry-independent kernels for computing the 
interior vector potential and its derivatives. All the ad- 
vantages demonstrated in this paper for surface methods 
will be retained. 



y 




FIG. 26: (Color) Illustration of a bent box with straight ends, 
used for computing a transfer map for bending dipoles with 
large sagitta. A design orbit is illustrated in red. 
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APPENDIX 

MaryLie/IMPACT is a hybrid code that utihzes 
Lie-algebraic niethods for computing and manipulating 
charged-particle transfer maps through 5"^ order, while 
space-charge effects are treated using Particle-in-Cell 
methods [11]. We made use of its Lie algebraic capa- 
bilities. In the Lie algebraic approach, maps are com- 
puted and manipulated in Lie algebraic form. Each map 
describes the transformation of the full six-dimensional 
phase-space coordinates of a particle as it passes through 
a given beam-line element. Because of the symplectic na- 
ture of Hamiltonian motion, through aberrations of order 
(n — 1) such a map has the Lie representation 

M=n2 exp(: h ■■) exp(: h :) ■ ■ ■ exp(: /„ :) (62) 

where TZ2 describes the linear part of the map, and each 
fj is a homogeneous polynomial of degree j describing 
the nonlinear part of the map [H [T^ . 

The linear map 1^2 and the Lie generators fj are de- 
termined by solving the equation of motion 



where 



M=M:-H 



H = H2 + H3 + Hi + -' 



(63) 



(64) 



is the charged-particle Hamiltonian expressed in terms of 
deviation variables about the design orbit and expanded 
in a homogeneous polynomial series. The deviation vari- 
able Hamiltonian H is determined in turn by the Hamil- 
tonian K for which path length is the independent vari- 
able. In Cartesian coordinates and with z taken as the 
independent variable, K is given by the relation 

K = - [{pt + q'^f/c^ - rr?c^ - (p, - qA^f 

- {py - qAy)^]'/^ - qA,. (65) 

Here $ and A are the electric scalar and magnetic vector 
potentials, respectively. 

We conclude that (in the case of no electric fields, 
$ = 0) what is needed are Taylor expansions for the 



components of A in the deviation variables x and y. In 
the straight-element case these Taylor expansions can be 
found using the circular cylinder harmonic based expan- 
sion ( 12 1. Suppose, for example, we wish to retain in the 



expansion of the Hamiltonian H appearing in ( 64 ) homo- 
geneous polynomials through degree 8. That is what is 
required to compute transfer maps through 7"^ order. If 
the design orbit lies on the z axis, as will be the case for 
any straight niultipole such as a solenoid, quadrupole, 
sextupole, octupole, etc., this expansion is straightfor- 
ward because in this case the Cartesian coordinates x, y 
are already deviation variables. We see from (651 that 



we must retain homogeneous polynomials in the vari- 
ables x,y through degree 7 in the expansions of A^ and 
Ay, and homogeneous polynomials in the variables x,y 
through degree 8 in the expansion of A^. Inspection of 
(12 1 shows that for the cases to = or to odd we then 
need the Cm]a{z) with (m + n) < 7. For the cases of even 
TO we need the Cm]a{z) with (to -|- n) < 8. 



In the case of a wiggler, the design orbit oscillates 
around the z axis. See Fig. 25. Now the requirement 



that the Hamiltonian H appearing in (64) be an expan- 



sion in deviation variables about the design orbit is more 
involved: The components of A must be expanded about 
the design orbit. As they stand in (12 1, they are ex- 
panded about the z axis. If we are to retain terms in H 
through degree 8 when a re-expansion is made about the 

In] 

design orbit, then coefficients Cm'a beyond those listed 
in the previous paragraph must in principle be included 
in the calculation. We say that higher-degree terms pro- 
duce lower-degree terms due to feed down. However, the 
effect of feed down is small if the oscillations are of small 
amplitude, as they are for the proposed ILC wiggler. In 
this example we have found that the approximation of re- 
taining only the terms in A^, Ay, and A^ through degree 
8 produces changes in the design orbit and the transfer 
map about the design orbit that are comparable to or 
smaller than the errors found in the benchmark studies 
of Section IV. We also remark that this feed-down prob- 
lem does not arise when the general geometry methods 
to be described in Part II are employed. 
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